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This  paper  presents  an  iterative  method  to  obtain  a  partial  or 
complete  solution  of  the  general  eigenproblem  which  does  not  require 
any  preliminary  modification  to  put  it  into  the  form  of  a  special  eigen¬ 
value  problem.  The  Rayleigh  Quotient  is  minimized  by  the  use  of  the 
conjugate  gradient  method  to  obtain  the  lowest  eigenvalue  and  the 
associated  eigenvector.  The  approach,  originally  proposed  by  Bradbury 
and  Fletcher  (Reference  1),  is  extended  to  permit  the  intermediate 
eigenvalues  and  eigenvectors  to  be  obtained  by  adapting  a  projection 
scheme  which  is  akin  to  Rosen’s  Gradient  Projection  Method  (Ref¬ 
erence  2).  This  technique  constrains  the  minimization  search  to  the 
subspace  M-orthogonal  to  the  previously  determined  eigenvectors.  A 
Theoretical  justification  is  presented  that  the  quadratic  convergence  of 
the  conjugate  gradient  method  is  preserved.  The  important  computer 
storage  advantages  of  the  conjugate  gradient  method  are  extended  by 
eliminating  the  need  for  assembled  stiffness  and  mass  matrices.  A 
number  of  structural  examples  are  presented  to  demonstrate  the 
effectiveness,  generality  and  stability  of  the  method. 
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SECTION  I 
INTRODUCTION 


A  frequent  intermediate  step  in  the  linear  dynamic  analysis  of  complex  structures  is  the 
solution  of  the  eigenproblem: 

KX  =  X  MX  (I) 

This  requires  the  determination  of  the  scalar  quantities  X  (eigenvalues)  and  the  corresponding 
non -trivial  solutions  X  (eigenvectors)  for  the  given  nxn  matrices  K  and  M  .  In  the  common 
structural  application,  K  and  M  are  respectively,  the  master  stiffness  and  mass  matrices 
of  the  structure,  and  their  order,  n,  corresponds  to  the  elastic  degrees  of  freedom  of  the 
system.  In  this  paper,  the  K  and  M  matrices  are  assumed  to  result  from  a  finite  element 
idealization  of  the  actual  structure. 

Frequently,  in  applying  the  formulation  of  Equation  1  to  the  study  of  the  vibration  character¬ 
istics  of  a  structure,  the  order  of  the  K  and  M  matrices  is  so  high  that  it  is  impractical  or 
prohibitively  expensive  to  obtain  the  complete  eigensolution.  On  the  other  hand,  to  carry  out 
a  reasonably  accurate  dynamic  analysis  of  the  structure,  it  is  possible  to  get  along  with  only 
a  partial  eigensolution.  It  is  this  class  of  problems  for  which  the  method  described  in  this 
paper  is  especially  useful. 

There  are  two  general  types  of  methods  for  the  eigensolution  of  Equation  1 :  transformation 
methods  and  iterative  methods.  The  transformation  methods  such  as  the  Jacobi,  Givens,  and 
Householder  schemes  (References);. are  almost  always  preferable  when  a  complete  eigen¬ 
solution  is  required.  On  the  other  hand,  the  labor  saving  involved  in  obtaining  only  a  partial 
solution  by  these  methods  can  be  a  small  fraction  of  the  total.  Furthermore,  the  transformation 
methods  accomplish  the  solution  by  operating  on  the  matrices  of  the  system  which  necessitates 
the  storage  of  large  matrices.  With  the  emergence  of  the  “consistent  mass  matrix,”  (Reference 
4)  another  difficulty  develops  because  of  the  necessity  to  transform  Equation  1,  the  general 
eigenproblem,  into  a  special  eigenproblem: 

DY  =  XY  (2) 

as  required  by  the  transformation  methods.  If  K  or  M  happen  to  be  sparse  or  banded,  this 
step  generally  produces  a  dynamic  matrix  D  with  more  extensive  storage  requirements  than 
either  K  or  M  . 
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The  direct  iterative  methods,  on  the  other  hand,  can  avoid  the  necessity  of  storing  the 
entire  matrix  by  using  modern  matrix  interpretive  methods,  yet  they  are  disadvantaged  in 
the  general  eigenproblem,  Equation  1,  and  generally  require  preliminary  modification  similar 
to  the  transformation  methods.  In  addition,  these  methods  are  plagued  by  convergence  dif¬ 
ficulties  and  are  computationally  expensive  for  the  intermediate  eigenvalues  and  eigenvectors. 


This  paper  describes  an  iterative  method  which  can  be  applied  directly  to  the  eigenproblem, 
Equation  1,  without  preliminary  modification.  It  uses  the  well  known  property  of  the  Rayleigh 
quotient. 


R  (  X  ) 


XT  KX 
XT  M  X 


(3) 


that  it  equals  the  eigenvalue  when  the  eigenvector  is  substituted  into  it  and  that  it  is  stationary 
in  the  neighborhood  of  an  eigenvector. 


The  basic  algorithm  is  simple:  The  Rayleigh  quotient  is  minimized  to  obtain  the  lowest 
eigenvalue  and  the  associated  eigenvector.  This  minimization  is  done  numerically  using  the 
conjugate  gradient  method.  Next,  for  the  second  eigenvalue  and  eigenvector,  the  Rayleigh 
quotient  is  again  minimized,  only  this  time  in  a  subspace  which  is  M-orthogonal  to  the  first 
eigenvector.  This  process  can  be  repeated  as  many  times  as  desired  to  obtain  as  many  of  the 
eigenvalues  and  eigenvectors  as  are  desired  up  to  the  complete  eigensolution.  This  approach  of 
obtaining  the  lowest  (or  highest)  eigenvalue  and  the  associated  eigenvector  was  originally 
proposed  by  Bradbury  and  Fletcher  (Reference  1). 

The  contribution  of  the  present  paper  is:  (1)  that  it  extends  the  approach  so  that  it  is 
practical  to  obtain  the  intermediate  eigenvalues  and  eigenvectors  without  a  lessening  of  the 
storage  and  efficiency  advantages,  and  (2)  it  explores  it  in  application  to  the  special  character¬ 
istics  of  finite  element  structural  dynamics  problems. 


The  extension  to  intermediate  eigenvalues  is  accomplished  by  using  a  gradient  projection 
scheme  (Reference  2)  for  constraining  the  minimization  search  to  the  subspace  M-orthogonal 
to  the  previously  determined  eigenvectors. 


The  advantages  in  structural  problems  of  the  formulation  of  Equation  3  are  that  both  the 
numerator  and  denominator,  as  well  as  all  of  the  other  quantities  required  by  the  iteration 
procedure  for  all  of  the  eigenvalues  desired,  can  be  computed  without  having  the  assembled 
K  and  M  matrices  at  hand.  This  is  accomplished  by  noting  that  the  numerator  is  twice  the 
strain  energy  for  a  given  X  (the  generalized  displacements),  and  that  the  denominator  is 
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twice  the  maximum  kinetic  energy  of  the  structure  and  that  these  can  be  computed  by  summing 
the  potential  and  kinetic  energies  of  the  individual  elements  of  the  discretized  structural  model. 
In  this  sense,  the  approach  is  an  extension  of  the  energy  search  method  documented  in  the 
literature  (References  5  and  6). 

Because  the  method  is  iterative  and  converges  quite  rapidly  when  reasonable  initial 
estimates  of  the  eigenvector  are  available,  it  lends  itself  well  to  embedment  within  structural 
optimization  procedures  where  dynamic  behavior  is  to  be  considered.  This  is  because  as  the 
optimum  structural  design  is  evolved,  the  eigensolution  generally  is  expected  to  change  only 
incrementally  from  design  to  design.  Hence,  the  previous  solution  provides  good  initial 
estimates  of  the  eigenvectors. 


SECTION  II 

FORMULATION  OF  THE  PROBLEM 

The  eigenproblem,  Equation  1,  can  be  written  as 

[k  -xm]  X  =  0  (4) 

If  X  is  a  solution  to  Equation  4  .then  bX  is  also  a  solution  for  any  nonzero  value  of  the  scalar 
b,  thus  the  eigenvector  corresponding  to  any  eigenvalue  X  is  arbitrary  to  the  extent  of  a  scalar 
multiplier.  In  other  words,  the  Rayleigh  quotient  defined  in  Equation  3  has  no  unique  minimum, 
but  takes  on  the  same  value  at  every  point  along  any  line  in  the  n-dimensional  space  passing 
through  the  origin.  Furthermore,  the  quotient  is  not  defined  at  the  origin.  Consequently,  the 
minimization  of  the  Rayleigh  quotient  is  not  quite  as  simple  as  that  of  a  function  with  a  well 
defined  minimum. 

The  redundant  degree  of  freedom,  which  prevents  us  from  determining  the  absolute 
magnitudes  of  the  components  of  the  eigenvector,  can  be  eliminated  by  an  arbitrary  normal¬ 
ization.  The  simplest  normalization  for  the  present  purpose  is  to  set  any  non-zero  component 
of  the  eigenvector  equal  to  one. 

The  Rayleigh  quotient.  Equation  3,  equals  the  eigenvalue  when  the  eigenvector  is  sub¬ 
stituted  in  it.  Moreover,  it  is  stationary  in  the  neighborhood  of  an  eigenvector  and  its  value  is 
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bounded  by  the  lowest  and  highest  eigenvalues  of  the  physical  system.  Thus,  the  minimization 
of  the  Rayleigh  quotient  will  yield  the  lowest  eigenvalue. 


The  minimization  problem  to  find  the  lowest  eigenvalue  can  thus  be  stated  as: 


Find  X  =  X,  such  that 


XT  K  X. 

R(X,)  =  — - 

X  M  X, 


is  mi  n  imu  m  ,  su  bject  to 


s  X,  «q  =1 


where  x  is  the  normalizing  or  reference  component  and  e  is  a  vector  with  its  qth  com- 
q  i  q 

ponent  as  one  and  zero  elsewhere  (i.e,  is  a  unit  coordinate  vector  for  the  qth  coordinate). 


An  illustrative  example  with  a  geometrical  interpretation  might  be  convenient  to  elucidate 
the  underlying  idea.  Consider  the  three  degrees  of  freedom  system  depicted  in  Figure  1. 
The  tubular  member  AB  of  mean  diameter  D  =  0,8  in.  and  the  wall  thickness  t  =  0.2  in.  is 
held  fixed  at  the  end  A  and  is  hinged  at  the  end  B.  A  model  consisting  of  two  standard  beam 
elements  was  used  and  the  vertical  displacement  of  the  middle  point  C  of  the  beam  and  the 
rotation  at  the  points  C  and  B  are  taken  to  be  the  three  degrees  of  freedom. 


Mean  diameter  (D)  =  0.8" 

Wall  thickness  (t)=  0.2" 

L  «  16.16" 

E=30  x  |06 1 bs./in? 
/>  =  0.28  lbs. /in3 

Figure  1.  Tubular  Beam  Fixed  at  the  End  A  and  Hinged  at  the  End  B. 
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The  stiffness  and  the  mass  matrices  of  the  structure  are  given  as: 
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Let  the  vector  space,  comprised  of  the  degrees  of  freedom  be  represented  as 


X  = 


*3 


and  if  we  pick  the  normalizing  component  as  x(  ,  then  the  normalized  vector  has  only  two 
unknown  components: 


X  s 


where  o  = 


x2 


and  b  - 


*3 


x2A, 

x3/x. 


I 

a 

b 


Note  that  in  this  example  we  have  chosen  e  ^  vector  as 


I 

0 

0 


From  the  expressions  of  K  ,  M  and  X  given,  we  can  write  the  Rayleigh  quotient,  R,  ex¬ 
plicitly  as: 

2  . 


0  73~6.0b  +  63.5a  +  5l.7ob  4-  31.7b  x  |Q 
0.44  4- 0.6b  4- 2.9a2  -  2. 2ab  +  1.5  b2 


A  plot  of  R(  X  )  for  different  values  of  a  and  b  is  shown  in  Figure  2,  which  represents  the 
contour  map  of  the  values  of  the  Rayleigh  quotient  corresponding  to  the  normalized  modes  of 
the  system.  As  is  seen  in  Figure  2,  the  Rayleigh  quotient  takes  on  the  minimum  value  at  the 
point  1,  a  maximum  value  at  the  point  3  and  an  intermediate  value  at  the  saddle  point  2. 
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Figure  2.  Contour  Map  of  Rayleigh  Quotient. 
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The  first  three  eigenvalues  X(  ,  \2  ,  and  X3  and  their  associated  eigenvectors  X,  . 
X  2  and  X  3  are  given  by 

X,  i  0.780X10®  ,  X=  0.  I099X  10®  ,  X=  0.780  X 10®  (rodions  /sec.)2 

I  Z  o 

and 


1.0  \ 
-0.033 

0.122  I 


I .0  \  /  1.0 

0.592  J  -  X3:  (  -0-793 

-0  662  /  \  -  I  .29 


The  minimization  problem  posed  in  Equations  5,  6  will  yeild  X  and  X.  . 


In  problems  of  structural  dynamics,  the  eigenvectors  represent  the  mode  shapes  and  the 
choice  of  a  non-zero  component  ordinarily  presents  no  serious  problem.  The  mechanization 
of  this  aspect  of  the  method  is  briefly  described  in  the  discussion  of  Example  1,  Section  VII 
Numerical  Examples. 


Once  the  lowest  eigenvalue  X(  =R(  X( )  is  known,  the  next  higher  or  (second)  eigenvalue 
and  its  associated  eigenvector  can  be  determined  by  posing  a  new  minimization  problem: 


Find  X  =  X2  such  that 


is  minimum,  subject  to 


and 


R{X2)  : 


x2  kx2 

XIM*2 


Xl*r 


X2  M  X,  =  0 


(7) 

(8) 

(9) 


In  the  subspace  defined  by  the  constraints.  Equations  8  and  9,  the  Rayleigh  quotient  takes 
on  a  unique  minimum  (assuming  distinct  eigenvalues)  at  the  eigenvector  associated  with  the 
second  lowest  eigenvalue.  The  constraint  Equation  8  is  of  the  type  already  discussed  and 
9  represents  the  imposition  of  the  M-orthogonality  condition  between  X,  and  X  2.  Geometrically 
speaking,  these  constraints  merely  restrict  the  portion  of  vector  space  in  which  the  search  for 
the  second  eigenvector  is  carried  out  and  in  this  restricted  space  R  has  a  minimum  cor¬ 
responding  to  X 
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The  determination  of  the  third  and  subsequent  eigenvalues  together  with  their  associated 
eigenvectors  up  to  the  complete  eigensolution  is  accomplished  by  solving  a  sequence  of  problems 
similar  to  the  one  presented  by  Equations  7  to  8.  The  only  change  is  that  each  time  one 
additional  equation  of  constraint  has  to  be  imposed  on  the  minimization  problem  to  satisfy  the 
condition  that  the  eigenvector  currently  being  sought  is  M-orthogonal  to  all  of  the  previously 
determined  eigenvectors.  The  problem  of  determining  the  ^th  eigenvalue  (2  <  £  <  n)  can  thus 
be  written 


Find  X  =  Xi  such  that 


R(X£1  = 


KXi 

xiMXi 


subject  to 


H'i 


-  I 


(10) 


(II) 


and 

Xj  M  X.  =  0  ,  i  =  1,2, ■■  ■',/  =1  (12) 

where  X.,  i  =  1,2,..., £-1  are  assumed  to  be  known  when  the  Jlth.  eigenvector  is  being 
sought. 


Denoting 


MX  =  V 

i  i 

the  constraint  Equation  12  can  be  written  as 

x\  Vj  --  0,  i  =  I  ,  2  ,  •  •  ,  X-  \ 


(13) 


(14) 
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SECTION  III 

MINIMIZATION  ALGORITHM 


The  methods  of  minimizing  a  function  of  several  variables  can,  in  general,  be  classified 
as  gradient  methods  and  non-gradient  methods.  The  gradient  methods  use  the  local  information 
about  the  rate  of  change  of  the  function  with  respect  to  the  changes  in  the  variables  and  require 
the  evaluation  of  the  gradient  vector,  inthiscase  VR.  These  methods  are  inherently  the  more 
powerful  as  more  information  about  the  function  is  used  and  are  preferred  over  the  nongradient 
methods. 


The  Rayleigh  quotient  as  a  function  of  the  n  variables  (x, ,  x2,... ,xn)  =  X  defined  in 
Equation  3  is  differentiable,  and  its  gradient  vector 


VR  =  g  = 


2  K  X 
XT  MX 


(  XT  K  X  ) 


(XT  MX) 


2(  KX  -  R  M  X) 
(XT  M  X  ) 


(15) 


is  easily  computed.  Therefore,  it  is  logical  to  carry  out  the  minimization  of  the  Rayleigh 
quotient  by  one  of  the  gradient  methods. 


Through  the  middle  1950's,  one  of  the  most  popular  gradient  methods  was  the  method  of 
steepest  descent.  This  method  chooses  each  direction  of  search  to  be  the  negative  of  the 
gradient  vector.  Though  used  moderate  success  on  a  variety  of  problems,  it  often  turns  out 
to  be  hopelessly  slow  because  of  the  fact  that  successive  moves  are  perpendicular  to  each 
other  and  the  method  gradually  settles  into  a  steady  n-dimensional  zig-zag  for  functions  having 
any  significant  eccentricity.  The  convergence  difficulties  of  the  steepest  descent  method  have 
been  largely  eliminated  by  a  modification  of  the  basic  iteration  which  has  been  called  the 
conjugate  gradient  method  (Reference  7).  This  method  has  the  property  that,  for  a  quadratic 
function  of  n  variables,  it  will  converge  in  n  steps,  apart  from  round  off  errors.  For  general 
functions,  as  the  iterate  approaches  the  minimum,  the  function  is  usually  more  nearly 
approximated  by  a  quadratic  and  so  convergence  accelerates  toward  the  solution. 
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The  method  of  Davidon  (Reference  8)  (1959)  which  was  amended  by  Fletcher  and  Powell 
(Reference  9)  (1963)  is  also  a  gradient  technique  which  has  the  property  of  quadratic  con¬ 
vergence.  However,  the  simplicity  of  the  conjugate  gradient  method  and  its  modest  demands  on 
storage,  compounded  by  the  successful  experience  of  Fletcher  and  Bradbury  (Reference  1). 
dictates  the  use  of  conjugate  gradient  method  to  minimize  the  Rayleigh  quotient. 

The  algorithm  can  be  written  as : 

XQ  r  arbilrory  (a) 


Go 

=  VR(X0) 

(  b) 

S 

0 

'  "®o 

(  c  ) 

xi+. 

■  xi+a>s( 

(  d) 

(16) 

Gm 

*  7R,Xi<.|> 

(  e  ) 

/3j 

(  f  ) 

si+, 

(9) 

where  the  step  length  a  Ms  the  value  of  a  which  minimizes  R(  X  .  +aS.).  From  Equation 
16g  we  note  that  S  ,+1  is  alinear  combination  of  6  ,+1  and  Sq,  S(  S.  and  hence  it  is  a 

linear  combination  of  G  Q»  G,  G  i+r  The  algorithm  is  based  on  a  Gram-Schmidt  orthog- 

onalization  of  the  G.  and  the  derivation  is  documented  in  the  literature  (Reference  10). 

The  method  described  in  Equations  16  is  applicable  in  principle  to  any  unconstrained 
minimization  problem.  It  will  be  noted  that  the  constraint  Equation  6  is  trivially  satisfied  by 
setting  the  qth  component  of  the  starting  point  to  be  unity  and  setting  the  corresponding  com¬ 
ponent  of  the  gradient  vector  to  be  zero  throughout  the  search  space.  Thus,  the  problem  of 
minimizing  the  Rayleigh  quotient  function  to  find  the  lowest  eigenvalue  is  similar  to  an  un¬ 
constrained  minimization  problem  and  the  Fletcher  Reeves  algorithm  can  be  directly  applied. 

However,  the  use  of  conjugate  gradient  method  for  finding  the  intermediate  eigenvalue  is 
possible  only  when  the  minimization  of  the  Rayleigh  quotient  is  restricted  to  a  subspace  of 
X  in  which  the  constraint  Equations  11  and  12  are  continuously  satisfied.  In  order  to  insure  that 
the  search  is  carried  out  in  the  desired  subspace  of  X  ,  it  is  necessary:  (1)  to  start  the 
iteration  with  a  point  in  that  subspace;  and  (2)  to  project  the  gradient  vector  g  ,  Equation  15 
on  to  that  subspace.  Both  of  these  requirements  necessitate  the  use  of  some  matrices  which 
project  the  arbitrary  startingpoint  and  the  gradient  vector  g  on  to  the  subspace  of  constraints. 
Henceforth,  such  a  matrix  will  be  called  the  “projection”  matrix  and  the  way  it  is  generated 
is  discussed  in  the  subsequent  section. 
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SECTION  IV 


PROJECTION  MATRIX 


Let  P  be  a  matrix  which  has  the  property  that  for  any  vector  W  ,  the  vector 


satisfies 


Wp  =  P  W 


WpZj  =  0  ,  i=  I,  2,-  •  •  •  ,  q 


(  17) 

(18) 


where  Z  i  =  l,2,...,q  are  q  linearly  independent  vectors. 

Note  that  Equation  18  can  also  be  written  in  matrix  form  as: 

nt  Wp  =  0 


where 


(q  x  n  )  (n  xt ) 

(n  xq  ) 


(19) 


(20) 


In  other  words,  the  matrix  operator  P  eliminates  from  the  vector  W  the  non -orthogonal 
components,  thus  giving  the  vector  W  which  is  orthogonal  to  the  subspace  spanned  by  the 

Jt“* 

vectors  Z  i  =  1,2,. ..,q.  This  idea  can  thus  be  expressed  differently  as: 

1  q 


or,  in  the  matrix  form  as: 


Wp  =  W  -  I  u.  Z. 
i  =  I 


Wp  *  W -  N  U 

(n  xq)  (q  x  I  ) 


(21  ) 
(22) 


where  the  components  of  vector  U  are  u.,  i  =  1,2 . q.  Premultiplying  Equation  22  by 


N  we  obtain  from  Equation  19 


NT  Wp  s  NT  W  -  (  NT  N)U  =  0 


Therefore 


U  z  (  NT  N>"  ntw 
From  Equations  24  and  22  we  obtain 

Wp  =  W  -  N  (  NT  N )"'  NT  W 


(2  3) 


(24) 


=  {  I  -  N  (  NT  N)"'  NT}w 


(25) 


where  I  is  the  identity  matrix. 
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T 

Note  that  (  N  N  )  will  be  a  (qxq)  symmetric  matrix,  and  is  nonsingular  since  N  is  a 
(nxq)  matrix  composed  of  q  linearly  independent  columns.  Hence  (  NTN  )"’  exists. 

Comparing  Equations  17  and  25  we  obtain  the  projection  matrix  as: 

p  «  {  I  N  (  NT  N  )"'  NT  }  (26i 

In  particular,  for  determining  the  £  th  eigensolution,  a  projection  matrix 


pi  =  {  I  ~  *£  <  n}  Nj  )-'  Nj  }  (27) 

where  J 

=  [*j  ■ v,  -v y^-i  1  i28i 


will  project  the  gradient  vector  fl  ,  Equation  15  on  to  the  subspace  of  constraints  defined  by 
Equations  11  and  12.  Note  that  the  column  vectors  of  the  matrix  N £  are  linearly  independent. 
Thus 

9P  ^  P£  9  (29) 

A  question  which  immediately  warrants  attention  is,  since  we  use  a  instead  of  a  . 

♦  P 

will  the  conjugate  gradient  method  actually  converge?  In  other  words,  is  there  a  function  of 

which  0  p  is  the  gradient  and  whose  function  value  equals  R  in  the  subspace  defined  by 

Equations  11  and  12.  To  see  that  the  answer  is  affirmative,  it  is  only  necessary  to  construct 

the  Lagrange  function  for  the  constrained  minimization  problem,  Equations  10  to  12.  Define 

a  ‘ ‘ Lagrange-Rayleigh’ ’  function  for  the  ^th  eigensolution  as: 

Li  =  R(  Xi)_  ui(X/ej  -l)  “  2  UiX^V;_,  (30) 

i-z 

where  u.,  i  =  1,2,...,  X  are  the  Lagrange  Multipliers.  Note  that  the  stationary  condition  over 
the  variables  X  is 


l 

VL^VR-u,.  -  I  u,VM  =  0  ,3,1 

i  s  2 

or  in  matrix  form 

VLii  «  ~  Ul  U  =  0  (32) 

and  if  the  Lagrange  multiplier  vector  U  is  defined  by  Equation  24  as 

U  *  (  N^f*  nJ  o  (33) 
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we  obtain 

VLf  =  9  -  NjlUlJ  N^)'1  N]  g  s  g  =0p  (34) 

where  Pj)  is  given  by  Equation  27.  Furthermore,  L  =  R  in  the  proper  subspace,  because  the 
second  and  third  terms  of  Equation  30  are  identically  zero  by  virtue  of  constraint  Equations 
11  and  12. 

Thus,  the  “  Lagrange-Rayleigh”  function  has  the  same  value  as  the  Rayleigh  quotient  in 
the  proper  subspace  and  the  projection  matrix  P j  projects  the  gradient  vector  g  on  to  this 
subspace,  giving  thereby  the  gradient  of  the  “Lagrange-Rayleigh”  function. 


SECTION  V 

RECURSION  PROCEDURE 


It  is  always  possible  to  form  the  (£x£)  matrix  ( H  £  )  and  then  invert  it  to  obtain  the 
projection  matrix  P J  ,  Equation  27,  when  the  search  for  the  _^th  eigenvalue  and  its  cor¬ 
responding  eigenvector  is  made.  However,  it  is  desirable  to  avoid  this  computation,  as  the 
value  of  J?  will  be  large  when  higher  eigenvalues  are  searched. 


It  will  be  noticed  from  Equation  28  that  the  size  of  the  rectangular  matrix  N £  increases 
by  one  column,  every  time  a  new  eigenvalue  and  its  associated  eigenvector  is  searched.  This 
suggests  that  some  sort  of  a  recursion  procedure  should  be  used  which  permits  the  insertion 
of  the  vector  ,  on  to  the  set 


(35) 


and  uses  (  N  ,  Hj_  (  )''  .  which  is  presumed 
obtain  (  N Nj  )"'  from  (  N^_,  )  ' 


to  be  known.  Such  a  recursion  procedure  to 
is  described  in  Reference  2.  It  is  based  on 
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the  method  for  the  inverse  of  a  matrix  in  terms  of  the  inverses  of  certain  of  its  partitions. 
For  purposes  of  the  present  work,  let 


V."  1  Ni-, 


Therefore, 


Ni-,' 


eT  e. 

J  J 

•I*, 


.V 

J  2 


•I  y2 


SYM 


Vi  V, 


vT  V 

I  2 


VN-S 


vTv 

2  2 


h  -  1*JL  )  ■ 


vi-!v/-2 


— 

1 

!  *T  Yj 

1  j  U  —  \ 

1 

1 

!  v^- 

1 

1 

1-  _ 

> 

r 

>- 

■  j 

(36) 


where 


A,  =  Q 


We  presume  that  A' 


is  known,  then 


H  *  1  N]  Nx  1 


(37  ) 


(  38) 
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can  be  obtained  from 


and 


B 


B. 


1 


+  T<\ 


-i 


»  --f  (  A 


-i 


A  AT  A"1  1 

2  2  I  ' 


(39) 


8  3  ST 


where  the  scalar  quantity  s  is  given  by 


s  -  A ,  —  A„  A  . 


(40) 


A  procedure  to  obtain  B  ,  ,  B  2  ,  and  B  3  which  is  even  more  efficient  than  the  direct 
formulas  shown  above  and  which  uses  only  the  old  projection  matrix,  p£  >  the  old  inverse, 
Q  j'_  (  and  the  vector  V^_)  is  described  in  Appendix  I. 
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SECTION  VI 

STEP  BY  STEP  PROCEDURE 

CHOICE  OF  A  STARTING  POINT 

The  iterative  methods  should  have  a  good  starting  point,  otherwise  unnecessary  time  is 
wasted  inside  the  minimization  procedure  to  minimize  the  function.  However,  there  seems  to 
be  no  simple  method  for  finding  a  good  starting  value  for  any  iterative  method.  In  Reference  1, 
a  choice  of  the  starting  point  is  made  on  the  basis  that  the  first  step  of  the  minimization  pro¬ 
cedure  makes  the  fastest  descent  towards  the  minimum.  It  is  a  unit  vector  e  .  (the  jth  element 
of  e  .  is  S.j)  lying  along  the  coordinate  axes.  Such  a  starting  point  gives  a  badly  distorted 
mode  shape  and  was  not  found  to  be  the  best  choice  for  solving  the  structural  dynamics 
problem. 


Contrary  to  expectation,  a  set  of  random  vectors  proved  to  be  superior  to  the  selected 
unit  vectors  and  in  this  work  they  were  used  as  starting  points.  The  method  showed  reasonably 
good  convergence  from  these  points.  Once  the  solutions  converged,  a  knowledge  of  the  mode 
shapes  of  the  structural  system  was  obtained.  Experiments  were  conducted  in  which  some  of 
the  design  variables  of  the  structure  were  changed  and  the  eigensolution  of  the  modified 
system  was  obtained  by  using  the  mode  shapes  of  the  original  design  as  the  starting  points. 
This  showed  rapid  convergence. 


The  starting  point  for  the  search  of  the  lowest  eigenvalue  needs  to  satisfy  only  one 
constraint.  Equation  6,  which  is  trivially  satisfied  by  dividing  through  by  the  qth  component. 
However,  the  minimization  algorithm.  Equations  16,  generates  a  sequence  of  vectors  which, 
in  the  limit,  tend  directionally  to  the  minimum  eigenvalue  on  the  search  space.  Thus  the  qth 
component  of  the  vector  s  so  generated  have  to  be  maintained  as  unity,  so  as  to  satisfy  the 
constraint,  Equation  6,  continuously  in  the  space.  This  is  achieved  by  setting  the  qth  com¬ 
ponent  of  the  gradient  vector  at  the  particular  point  equal  to  zero.  In  other  words,  no  “move” 
is  made  in  the  qth  direction  of  the  search  space. 


The  starting  point  for  the  search  of  second  eigenvalue  has  to  satisfy  an  additional  con¬ 
straint,  Equation  9,  and  this  can  be  easily  satisfied  by  Schmidt  orthogonalization. 


Let 


be  some  initial  estimate.  Therefore, 


y(0)  .  tO) 

A  2  A  2 


V  )V 

¥l  ¥l 


(41  ) 
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satisfies  Equation  9  if  v}  V,  =  1.  Note  that  V(  *  MX,  ,  Equation  13  where  X!  is  the 
eigenvector  corresponding  to  the  first  (lowest)  eigenvalue. 

In  order  to  get  the  starting  point  for  the  search  of  subsequent  eigenvalues,  the  initial 
estimate  is  passed  through  a  projection  matrix  as  discussed  in  Section  IV,  Projection  Matrix, 
where  the  matrix  is  given  by  Equation  27  but  the  matrix  N  £  where 

"i  (421 

is  to  be  used  instead  of  matrix  N £  ,  Equation  28. 

FUNCTION  EVALUATION 

The  value  of  the  function  to  be  minimized  is  required  at  the  end  of  each  iteration  cycle 
for  almost  all  the  iteration  methods  since  the  convergence  criterion  is  based  on  the  function 
value.  Moreover,  for  the  particular  problem  of  minimizing  the  Rayleigh  quotient,  the  function 
value  is  required  at  each  cycle  in  order  to  evaluate  the  gradient  vector.  Equation  15.  It  is, 
therefore,  necessary  to  have  an  efficient  routine  for  function  evaluation,  in  order  to  avoid  the 
time  which  would  be  otherwise  wasted  inside  the  minimization  procedure.  The  value  of  the 
Rayleigh  quotient  can  be  computed  without  having  the  assembled  K  and  M  matrices  at  hand. 
This  is  accomplished  by  noting  that  the  numerator  is  twice  the  strain  energy  for  a  given  X  (the 
generalized  displacements) ,  and  that  the  denominator  is  twice  the  maximum  kinetic  energy 
of  the  structure  and  these  can  be  computed  by  summing  the  potential  and  kinetic  energies  of 
the  individual  elements  of  the  discretized  structural  model.  Thus  R(  X )  can  be  written  as 

y  yV  y. 

;-i  1  1  1 

R(X)  -  (43) 

r  J 

T  Y.  m  Y. 

^  i  i  i 

i  =  l 

where  r  is  the  number  of  discrete  elements,  k  i  and  m  .  are,  respectively,  the  stiffness  and 
mass  matrices  of  the  ith  element  of  the  discretized  structure  and  Y  .  is  the  displacement 
vector  of  the  ith  element  corresponding  to  the  generalized  displacement  vector  X  .  To  obtain 

the  vector  Y^,  i  =  1,2 . r  from  the  vector  X  is  rather  easy  and  is  a  matter  of  logical 

operations.  A  variety  of  methods  exist  for  such  logical  operations.  One  such  scheme  is 
described  in  Reference  11  and  an  equivalent  technique  is  discussed  in  detail  in  Appendix  II. 
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For  the  purposes  of  the  present  work,  a  separate  subroutine  was  written  in  Fortran  IV 
to  decompose  the  vector  X  to  the  vectors  Y  .,  to  evaluate  the  vectors  k.Yj  ajid 
m  .  Y.  and  then  to  assemble  back  the  product  KX=A  and  MX  5  8.  The  assembly  of  the  vectors 
A  and  B  is  accomplished  through  a  logical  operation  on  the  vectors  k.Yj  and  rrr  Y  . 
which  is  merely  the  inverse  of  the  operation  described  in  Appendix  II. 

As  will  be  seen  later,  this  subroutine  was  used  over  the  over  again  to  evaluate  certain 

other  quantities  required  in  the  iteration  process,  other  than  the  function  evaluation  which  is 

T  T 

now  obtained  by  two  vector  multiplications  X  A  and  X  B  and  one  division  in  order  to 
get: 

X^A 

R  (  X  )  = 

x'b 

Needless  to  emphasize  the  advantage  gained  by  getting  along  without  the  assembly  of  K  and 
M  matrices.  Frequently,  the  order  of  K  and  M  matrices  encountered  for  large  structures 
is  so  high  that  it  is  impractical  or  prohibitively  expensive  to  study  their  vibration  charac¬ 
teristics  without  making  approximations.  The  size  of  element  stiffness  and  mass  matrices 
are  relatively  much  smaller  than  the  size  of  the  assembled  stiffness  and  mass  matrices  of 
large  complex  structures.  Furthermore,  advantage  can  also  be  taken  of  similar  elements. 
For  example,  in  a  structural  system  of  a  large  number  of  elements,  if  only  three  types  of 
elements  are  used,  then  we  need  store  only  the  stiffness  and  mass  matrices  corresponding 
to  these  three  elements,  rather  than  for  all  of  the  r  elements. 

GRADIENT  EVALUATION 


The  minimization  algorithm,  Equation  16,  requires  the  evaluation  of  the  gradient  vector 
at  each  cycle  of  the  iteration.  The  gradient  vector,  g  ,  of  the  Rayleigh  quotient  function  is 
given  by  Equation  15.  For  determining  the  lowest  (first)  eigensolution,  the  component  of  the 
gradient  vector  corresponding  to  the  normalizing  component  of  the  eigenvector  is  set  equal  to 
zero,  in  order  to  continuously  satisfy  the  constraint  imposed  due  to  the  normalization  of 
the  eigenvector.  While  determining  the  intermediate  eigensolution,  the  projection  matrix 
*JL  ,  Equation  27  is  used  to  project  this  gradient  vector,  g  ,  on  to  the  proper  subspace  of 
search,  M-orthogonal  to  the  previously  determined  eigenvectors.  The  component  of  the  gradient 


vector,  g  ,  corresponding  to  the  normalizing  component  of  the  eigenvector  turns  out  to  be 
zero  automatically,  but  a  small  number  often  appears  due  to  roundoff  errors  and  this  is  simply 


removed  by  setting  that  component  equal  to  zero. 
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EVALUATION  OF  STEP  LENGTH 


Once  a  direction  of  move  S  .  has  been  chosen,  we  must  determine  a  *  so  that  the 
function  is  minimized  in  that  direction.  Thus  the  problem  of  finding  the  step  length  is  essentially 
the  linear  search  problem  which  requires  the  determination  of  the  ct  *  along  S  ^  through 
X  .  at  which  the  value  of  the  Rayleigh  quotient  function 

(X  fa.S  )T  K  (X  +  a  S  ) 

R  (  X .  +  a. S.  )  =  - 1 - -  (44) 

1  1  '  (X,  +  ctj S j)  M(X.  +  aiS] ) 

is  a  minimum,  i.e., 

5R((a*  )  =  0  (45) 

=  a*-  1 

i  i 

In  the  general  problem,  no  expression  is  available  to  determine  a  so  an  interpolation 
approach  (References  7,  8,  and  9)  is  adopted.  In  the  particular  problem  of  Rayleigh  quotient, 

however,  an  explicit  expression  in  can  be  generated  (Reference  1)  from  Equations  44 

and  45  which  has  the  form 

uai  +vaj  +  *  10 

where 

u  =  (ST  K  S.)(X[  M  S.)  —  (Xj  K  S.)(sj  MS.  ) 

v  =  (S[  K  Sj)(x{  M  Xj  )  -(  xj  K  Xj)(sT  M  S.  )  (47) 

w  =  (XT  K  Sj  )(xT  MX.  )  -  (x[  K  Xj  )(xT  MS.) 
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The  two  roots  of  Equation  46  correspond  to  the  maximal  and  minimal  points  of  the 
Rayleigh  quotient  in  the  direction  S.  through  X.  as  shown  in  Figure  3.  The  minimal 
function  value  corresponds  to  a?. 


R 


Note  that  the  matrix  products  required  in  Equation  47  to  evaluate  the  coefficients  of  the 
quadratic,  Equation  46,  are  easily  obtained  through  the  subroutine  discussed  above  in  sub¬ 
section  “Functional  Evaluation.” 
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SECTION  VII 

NUMERICAL  EXAMPLES 


In  all  the  illustrative  examples  presented  in  this  section,  (1)  a  general  planar  beam  element 
having  six  degrees  of  freedom  was  used  to  model  the  various  structures,  (2)  the  distributed 
mass  of  the  system  was  used  to  evaluate  the  components  of  the  mass  matrix  (the  “consistent” 
mass  matrix)  (Reference  4),  (3)  the  Univac  1108  digital  computer  was  used  to  obtain  the 
numerical  results  and  (4)  except  where  noted,  random  vectors  were  taken  as  starting  points 
for  the  minimization  algorithm  to  obtain  the  eigensolution. 

EXAMPLE  1 

As  a  simple  application,  the  cylindrical  cantilever  rod  shown  in  Figure  4a  is  considered. 
An  attempt  was  made  to  obtain  the  complete  eigensolution  of  this  simple  structure  by  the 
successive  minimization  of  the  Rayleigh  quotient  in  the  appropriate  subspace. 

Based  on  the  surmise  that  the  lowest  eigenvalue  would  correspond  to  the  first  cantilever 
mode,  the  degree  of  freedom  associated  with  the  transverse  deflection  of  the  cantilever  at  the 
free  and  (xQ)  was  made  to  equal  one  and  the  others  were  taken  as  zero  in  order  to  start  the 

O 

search  for  the  lowest  eigenvalue.  (This  obviously  makes  xg  the  normalizing  component).  A 
capability  was  built  into  the  computer  program  to  change  the  normalizing  component  whenever 
the  magnitude  of  any  other  component  of  the  vector  in  the  appropriate  subspace  exceeded 
five  times  the  magnitude  of  the  current  normalizing  component  (which  in  any  case  is  one). 
Every  change  of  the  normalizing  component  necessitated  the  restart  of  the  minimization 
algorithm  due  to  the  change  in  the  subspace  of  search. 

Curiously  enough  the  complete  eigensolution  so  obtained  gave  the  first  six  eigenvalues 

and  their  associated  eigenvectors  and  did  not  pick  up  the  remaining  three  eigenvalues.  This 

is  because  the  three  axial  degrees  of  freedom  are  uncoupled  from  the  six  translational  and 

rotational  degrees  of  freedom  of  the  cantilever  beam  and  with  the  particular  choice  of  starting 

point,  the  method  could  not  enter  into  the  subspace  in  which  the  eigenvalues  associated  with 

the  axial  modes  of  vibration  are  located.  If  instead  of  e  ,  the  starting  points  are  taken  as 

o 

e7>  only  the  eigenvalues  associated  with  the  three  axial  modes  of  vibration  are  obtained. 
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Furthermore,  if  the  same  structure  is  oriented  differently  as  shown  in  Figure  4b  where 
the  translational,  rotational  and  axial  degrees  of  freedom  are  all  coupled,  the  complete 
eigen  solution,  i.e.  all  the  nine  eigenvalues  and  the  corresponding  eigenvectors,  were  obtained 
from  the  starting  point  of  e  . 

O 

8  9  9 

The  first  three  eigenvalues  of  the  model  are:  0.2466  x  10  ,  0.1161  x  10  and  0.9747  x  10 

4  5  5 

which  correspond  to  the  frequencies  of  0.4966x10  ,  0.1079  x  10  and  0.3122  x  10  radians/sec. 

4  5  5 

These  compare  favorably  with  the  exact  values  of  0.4945  x  10  ,  0.1065  x  10  and  0.3100  x  10 
radians/sec.,  and  respectively  correspond  to  the  first  cantilever  beam  mode,  the  first  axial 
mode  and  the  second  cantilever  beam  mode. 

EXAMPLE  2 


Consider  as  another  illustrative  example  the  planar  frame  shown  in  Figure  5a  consisting 
of  tubular  members  pinned  together  at  the  nodes.  Each  member  was  modeled  with  a  single 
beam  element.  Thus ,  the  structure  has  fourteen  degrees  of  freedom.  A  complete  eigensolution 
(i.e.  all  the  14  eigenvalues  and  eigenvectors)  of  the  system  was  obtained  by  successive 
minimization  of  the  Rayleigh  quotient.  The  time  taken  was  4.17  seconds. 

EXAMPLE  3 


The  partial  eigensolution  of  a  56  degrees  of  freedom  system  shown  in  Figure  5b  was 
obtained  by  minimizing  the  Rayleigh  quotient.  Each  member  of  the  four  bay  planar  frame, 
pinned  at  the  nodes,  was  modeled  with  a  single  beam  element.  The  determination  of  ten  eigen¬ 
values  and  eigenvectors  took  77.91  seconds. 
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ii 


H-H 

TYP. 


Mean  Dia.  (D)  =  1.6" 
Wall  Thickness  (t)  =  O.l" 

E  =  30  x  I06  lbs. /in2 

P  -  0.28  lbs. /in2 


(b) 


Figure  5.  Planar  Frames:  (a)  Fourteen  Degrees  of  Freedom: 
(b)  Fifty- six  Degrees  of  Freedom. 


295 


AFFDL-TR-68-150 


296 


Figure  6.  Planar  Frame  (134  Degrees  of  Freedom,  If  Each  Member  is  Modeled 
of  One  Element). 
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EXAMPLE  4 

A  frame  structure  having  134  degrees  of  freedom  shown  in  Figure  6  was  analyzed  to  obtain 
the  first  five  eigenvalues  and  the  associated  eigenvectors.  Each  member  was  pinned  at  its 
nodes  and  modeled  with  only  one  element.  The  number  of  iterations  required  to  converge  to 
the  solution  up  to  the  eighth  decimal  place,  the  time  taken  to  obtain  each  eigenvalue  and  the 
numerical  result  obtained  are  given  in  Table  I.  The  eigenvectors  so  obtained  were  used  as  the 
starting  points  for  the  minimization  algorithm  to  obtain  the  eigensolution  of  a  changed  design 
of  the  same  structure.  The  change  in  the  original  design  was  brought  about  by  increasing  the 
mean  diameter  of  members  1  to  20  by  10%  and  by  decreasing  the  mean  diameter  of  members 
21  to  29  by  10%.  The  results  for  the  changed  design  are  also  given  in  Table  I.  The  total  time 
taken  to  obtain  the  partial  (first  five)  eigensolution  of  the  original  design  was  46.95  seconds 
while  that  for  the  changed  design  was  20.30  seconds.  Thus,  the  method  converges  rapidly  for 
the  changed  design.  The  reason  is  that  the  eigenvectors  of  the  original  design  provide  the 
reasonable  initial  estimates  of  the  eigenvectors  for  the  changed  design  and  thus  are  good 
starting  points  for  the  minimization  algorithm. 


It  would  be  observed  from  Table  I  that  the  first  two  eigenvalues  are  well-separated  while 
the  next  three  are  closely  spaced.  A  study  of  the  associated  eigenvectors  explains  this  behavior. 
The  first  two  eigenvalues,  respectively,  correspond  to  the  mode  shapes  schematically 
represented  in  Figures  7a  and  b.  The  next  three  closely  spaced  eigenvalues  correspond  to  the 
mode  shapes  where  the  overall  structure  frequency  is  dominated  by  the  individual  frequency 
of  the  long  members  of  the  structure  marked  34, 35, 44,  and  45.  This  interaction  of  the  individual 
member  frequency  with  the  overall  structure  frequency  is  responsible  for  the  three  closely 
spaced  eigenvalues. 

TABLE  I 

FIRST  FIVE  EIGENVALUES  OF  THE  ORIGINAL  AND  THE  CHANGED  DESIGN.  EIGENVALUES 


REQUIRED  TO  CONVERGE  UP  TO  THE  EIGHTH  DECIMAL  PLACE 


mm 

Original  Design 

Changed  Design 

o 

Iterations 

Time 

Value 

Iterations 

Time 

Value 

wEm 

(secs.) 

(secs .) 

i 

59 

5.3 

.62309556xl06 

26 

2.22 

.64067844xl06 

2 

56 

4.27 

.58762266x10? 

40 

3.04 

.60087324x107 

3 

205 

19.53 

.75223340x10? 

73 

6.11 

.75745767x10? 

4 

100 

8.36 

.77330454x10? 

46 

4.50 

.75729213x10? 

5 

105 

_ 

9.49 

.79875840x10? 

43 

4.43 

.75753798x10? 

298 


AFFDL-TR-68-150 


The  same  example  was  rerun  with  a  relaxed  convergence  criterion  {this  time  the  solution 
was  required  to  converge  to  only  the  sixth  decimal  place)  and  the  results  obtained  are  given 
in  Table  II.  The  total  time  taken  for  the  partial  eigensolution  of  the  original  design  is  24.34 
seconds  while  that  for  the  changed  design  is  7.97  seconds.  Thus,  there  is  a  considerable  saving 
of  time  in  choosing  a  less  stringent  convergence  criterion.  However,  although  the  error  does 
not  seem  to  propagate  badly  from  eigenvalue  to  eigenvalue,  the  eigenvectors  so  obtained  are 
not  as  accurate  as  those  obtained  by  the  more  stringent  convergence  criterion. 

TABLE  XI 


FIRST  FIVE  EIGENVALUES  OF  THE  ORIGINAL  AND  THE  CHANGED  DESIGN.  EIGENVALUES 
_ REQUIRED  TO  CONVERGE  UP  TO  THE  SIXTH  DECIMAL  PLACE 


No.  of 
Eigne- 

value 

Original  Design 

Changed  Design 

Iterations 

Time 

(secs.) 

Value 

Iterations 

Time 

(secs.) 

Value 

1 

45 

3.49 

.623109xl06 

6 

0.56 

. 640686xl06 

2 

45 

3.51 

.587632xl07 

9 

1.14 

. 600891xl07 

3 

50 

4.15 

.  752796xl07 

22 

2.13 

. 757763xl07 

4 

55 

5.55 

. 773182xl07 

27 

2.11 

. 757269xl07 

5 

76 

7.64 

.  798727xl07 

28 

2.03 

.757271xl07 

EXAMPLE  5 


When  each  of  its  members  is  modeled  by  two  beam  elements  the  number  of  degrees 
of  freedom  of  the  system  shown  in  Figure  6  is  281.  The  results  for  the  first  five  eigenvalues 
of  this  refined  model  are  given  in: 


TABLE  III 

RESULTS  OF  281  DEGREE-OF-FREEDOM  MODEL  OF  FIGURE  6. 


No. 

Iterations 

Time 
(secs . ) 

Value 

1 

417 

71.78 

.6207  x  106 

2 

301 

55.57 

.5291  x  107 

3 

614 

117.22 

.6295  x  107 

4 

222 

42.48 

.6371  x  107 

5 

297 

59.12 

.6643  x  107 
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As  anticipated,  the  refined  modeling  gave  a  better  correlation  between  the  closely  spaced 

eigenvalues  of  the  overall  structure  and  the  eigenvalue  corresponding  to  the  lowest  natural 

frequency  of  the  individual  members  having  a  length  of  21.21  in.  (any  of  the  ones  marked  34, 

35,  44  and  45  on  Figure  6).  The  lowest  natural  frequency  of  the  simply  supported  tubular 

4 

beam  of  mean  diameter  1.6  in.,  wall  thickness  0.1  in.  and  a  length  of  21.21  in. is  0.251  x  10 

7 

radians/sec.,  which  correspond  to  an  eigenvalue  of  0.6300  x  10  .  Thus,  the  closely  spaced 
eigenvalues  are  the  result  of  the  interaction  between  the  individual  member  frequency  and  the 
overall  structure  frequency. 


SECTION  vm 
CONCLUSION 


An  iterative  method  to  obtain  the  partial/completeeigensolution  of  a  general  eigenproblem 
arising  in  structural  dynamics  is  described.  The  method  does  not  require  preliminary  modifi¬ 
cation  to  put  the  general  eigenproblem  into  any  special  form.  It  has  been  found  to  be  effective 
for  the  partial  eigensolution  of  complex  structures  and  thus  is  useful  for  the  dynamic  analysis 
of  complex  structures. 

Since  the  method  converges  rapidly  when  reasonable  initial  estimates  of  the  eigenvector 
are  available,  it  lends  itself  well  to  embedment  within  structural  optimization  procedures  where 
dynamic  behavior  is  to  be  considered. 
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APPENDIX  I 


PROCEDURE  TO  OBTAIN  Qj'  FROM  Q/_\  " 


As  discussed  in  the  text  of  the  paper,  Q~J  is  given  by: 

I 


Q/'  =  (  Hj>  M,f 


where 


B, 


®  i  ;  b2 

K  |  _ 

"•  4.  —  ( 1  £  fj  4~l  ) 

1  ^  S  v  Al  H2  M2  M I 


b  =  - -1-  (a;1  a  ) 
2  s  I  2 


B,  =  — 


and 


Furthermore, 


and 


It  is  presumed  that 


is  known. 


s  ”  A  j  A  2  A ,  Aj 


A,  =  0^.,  =  Ng 


A,  .  .  ,Vje.  , 

-  0}'.,  »lNTi_,N^|  )- 


(48) 


(49) 


(50) 


C5t  ) 


(52) 


We  can  write 


where 


A 2  *'  Ni-t  Vi- . 


N 


'i- 


=  [*)  '  V.  'V2'-'Vi-] 


(53) 
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Therefore,  from  Equations  48,  49,  50,  and  51,  we  obtain 

s  =vLvi-, -vI-, (nL 


s  V ,  pi-,  v, 


where 


p^-,  -  {  1  -  ■/.,  N N/-, } 
is  the  “projection”  matrix. 


(54) 


(55) 


The  projection  matrix,  Equation  55,  has  a  neat  property  that  the  vectors  Pj  (  w  and 
are  orthogonal.  This  can  be  seen  by  mere  multiplication  of  the  two  vectors: 


1 

tH 

}  w  =  WT  { 

1  ~»£ 

-  ‘"j-, 

1 

1 

^4 

{ 

-  ,  '  "J. 

,  NiV" 

1 

Z 

WT  { 

< 

Ni- 

.>“  Ni- 

Ni- 

,  <N/-, 

na.  >" 

* 

(56) 


since  (  Hj_{  Hj  )'•  (  n}_(  N^JbX 


The  term  within  the  brackets  j  j  on  the  right  hand  side  of  Equation  56  is  identically 
zero  and  hence  we  obtain 


"T  pi-,  {T  -  p/-,}  w  =0 


Equation  57  can  be  rewritten  as 


WT  W  :  WT  P.  P  a  W 

X-\  X-\  1  -i 


P  W 

£-\ 


Therefore,  we  obtain  from  Equation  54 


s  5  1  pi-  V- 


(57) 


(58) 


(59) 


Note  that  the  projection  matrix  Pj_{  is  already  known  and  thus  to  obtain  s  from  Equation 
59  is  computationally  more  efficient  than  that  from  Equation  50. 
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If  we  denote 


*7  v(NV.  N/-ir'  NJ->  v^-'  l60> 

then  from  Equations  55  and  60  we  obtain 

V.  vi-,  =  vi-,  -  V.  '/->  160 

Furthermore,  from  Equations  52,  60,  61,  and  49  we  obtain 

B,  ■  I  »/-,  +T  ri-i  rJ-l 

B2  =  -Tri-i  1621 

and 

B3  =  T- 

Thus  the  procedure  to  obtain  Q  J  from  Q^_  t  can  be  summarized  as: 

1,  Compute  the  vector  f£_{  from  Equation  60  as  two  matrix  vector  multiplications. 

2.  Compute  the  scalar,  s,  from  Equation  61  and  59. 

3,  Compute  the  matrix  ,  the  vector  B2  and  the  scalar  B3  from  Equation  62. 

4.  Form  the  desired  Q^1  from  Equation  48 
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APPENDIX  II 

TRANSFORMATION  BETWEEN  REFERENCE  AND  LOCAL  COORDINATES 


Let  SR  be  an  rx(t+l)  table  whose  first  column  contains  the  element  numbers  arranged 
from  1  to  r  and  whose  rows  contain  in  the  remaining  t  columns  the  corresponding  local  con¬ 
sistent  coordinates  of  the  elements,  replaced  one  to  one  by  the  number  of  the  corresponding 
degree  of  freedom  of  the  structure  in  the  reference  coordinate  system.  The  ith  row  of  the 

SR  table  is  used  to  obtain  the  (txl)  vector  Y  ^  i  =  1,2 . r  from  the  (nxl)  vector  X  .  Y  i  is 

the  displacement  vector  ofthe  ith  element  corresponding  to  the  generalized  displacement  vector 
X  of  the  system.  As  an  illustrative  example,  consider  the  simple  cantilever  beam  shown  in 
Figure  4,  which  is  modeled  with  three  general  planar  beam  elements.  Thus,  the  number  of 
elements,  r,  is  3,  the  local  degrees  of  freedom  of  an  element  in  the  reference  coordinate  sys¬ 
tem,  t,  is  6,  and  the  number  of  degrees  of  freedom  of  the  system,  n,  is  9.  The  SR  table  for  this 
system  is  shown  below: 


Element 

No. 

Degrees  of 

Freedom 

Axial 

Displacement 

At  End 

Transverse 

Displacement 

At  End 

Rotation 

of  End 

P  Q 

P  Q 

P  Q 

1 

0  1 

0  2 

0  3 

2 

1  4 

2  5 

3  6 

3 

4  7 

5  8 

6  9 

Let  the  generalized  displacement  vector  be: 
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Thus  we  can  say  that 


i,e.,  the  vector  X  is  transformed  to  the  vector  Y  .  through  the  ith  row  of  the  SR  table,  where 
the  non-zero  components  of  the  vector  Y  are  the  components  of  the  vector  X  associated 
with  the  ith  element. 
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